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Abstract 

We propose a simple Monte Carlo method to calculate the interfacial free energy 7 between the 
substrate and the material. Using this method we investigate the interfacial free energys of the hard 
sphere fluid and solid phases near a smooth hard wall. According to the obtained interfacial free 
energys of the coexisting fluid and solid phases and the Young equation we are able to determine 
the contact angle with high accuracy, cos9 = 1.010(31), which indicates that a smooth hard wall 
can be wetted completely by the hard sphere crystal at the interface between the wall and the hard 
sphere fluid. 
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I. INTRODUCTION 



The interfacial free energy(IFE) 7 between different materials or different phases of a 
material plays crucial roles in many physical phenomena such as wetting and spreading 
An accurate result for 7 is necessary to understand these phenomena completely. The 
hard sphere system as a simple model system of condensed matter has been investigated 
extensively. By changing the density of the system, it undergoes a first order phase transfor- 
mation from the liquid phase to the FCC solid phase . An interesting problem of the hard 
sphere system related to the IFE is whether a smooth hard wall can be wetted completely 
by the hard sphere crystal in the density region of coexistence. The problem was studied by 
different groups in the past but a definite answer is still awaiting 



simulation 



4, l5j . Recent computer 



6| by Dijkstra provided a strong evidence of complete wetting. In that paper 
the author investigated the effect of wall separation on crystalline layers at the walls in a 
double-wall system and extracted the complete wetting information from the phenomenon 
of surface freezing (including complete wetting and capillary freezing). Alternatively, the 
direct way to study this wetting behavior is to calculate the cosine of the contact angle in 
terms of Young equation, i.e. cos 6* = (7^^ — 7stt,)/7/s, here the subscripts /, s and w de- 
note respectively the fluid, solid and wall, and the superscript c denotes the bulk coexisting 
phase. In the complete wetting state cos^ given by the above formula will be larger than 



unity. The fluid- wall IFE 7^^ can be estimated via density functional theory [7|, l8|, scaled 
)article theory Qj, molecular dynamic simulation lOj and Monte Carlo(MC) simulation 



However, it is unfeasible to calculate directly the solid- wall IFE 7^^^, using the 



above methods except for the MC simulation pj|. Recently, several groups have calculated 



;he crysta 

3,[i3,[i6i, 



/me 



17 



t IFE 7jg of a hard sphere system by different computer simulation methods 



18|. However, the accuracies of the contact angle 6 from t 



leir calculation are 



. Comparing with 
Q. Therefore, it 



too low to conclude definitely whether or not the wetting is complete [6 1 
7jg, the error of the contact angle 9 comes mainly from 7^^ and 7^^ [uJ, 
is desired to accurately calculate the IFEs, which provide us not only the confirming of the 
complete wetting scenario but also a benchmark to test other theoretical methods. 

In this paper we propose a simple Monte Carlo method of calculating the IFEs 7/^^, and 
'jsw, and obtain the IFEs of the coexisting hard sphere fluid and solid phase(i.e. 7^^ and 
7^^) with high accuracy, which nearly confirms the suggestion that at the smooth hard 



2 



wall-hard sphere fluid interface a complete wetting phase transition occurs very close to the 
bulk freezing transition. Our method is inspired by the cleaving method of Davidchack and 
Laird IJ], where they used two moving walls to transform the hard sphere system from 



the bulk phase into the phase confined completely by two walls and vice versa, and the 
IFE can be obtained by calculating the work done in the reversible process. However, as 
mentioned in Davidchack's work [14'], irreversibility is introduced during cleaving a fiuid 
system using constructed hard walls. The reason of irreversibility is that the system doesn't 
reach the equilibrium state in a typical run, which can be reduced by increasing the duration 
of equilibration run. The problem can be avoided by designing an auxiliary path which 
connects the confined and the bulk phase and the states along the path are all in equi 



state during the simulation. The goal is achieved 
and the so-called fiat histogram methods 0, 21 



ibrium 

22l |. From these method 'jf^ and jsw can 



combine the cleaving technique 



)e accurately estimated. In the following sections we employ Wang-Landau (WL) algorithm 

2l| to demonstrate the implementation of the current approach. 

In order to check the validity of our method, we calculated the IFE between the hard 
sphere system and the smooth hard wall at usual densities, including 7/„, and '-fsw for three 
different crystal orientations (1, 1, 1), (1, 1,0), and (l20,0). The results are consistent with 
those of Heni and Lowen within the statistical error |ll| . We also investigated the finite size 
effect of longitudinal dimension(along z direction) on the IFE, and then obtained the effective 
interaction potential between two hard walls. Finally, we use the method to calculate the 
IFEs of the coexisting phase and 7^^, and then obtained the contact angle cos6'. The 
paper is organized as follows: in section II, the model system and the numerical algorithm 
are described. The results are presented in section III. Section IV is a brief conclusion. 



II. MODEL AND ALGORITHM 



Our system is composed of hard sphere particles with diameter a in a fixed volume V. 
The periodic boundary conditions are imposed in all three directions. Hard sphere system is 
athermal, therefore its thermodynamic properties are solely dependent on the dimensionless 
packing fraction rj = vrpcr^/G, where p is the number density of the particles. Pairs of 
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particles i and j interact via the potential 







oo 



when I — I < 0", 
otherwise, 



(1) 



here r is the positions of the particles. 




ail a2 



FIG. 1: Diagram illustrating the change of the system boundary conditions by two synchronously 
moving smooth hard walls. When the walls are respectively placed on |ai| and \a2\ (|ai| = \a2 — 
Zq\ = Dq > cr/2), they do not interact with the particles and the system is in the bulk phase; when 
the walls are placed at the boundaries, the system are completely confined by two hard walls. 

Two moving smooth hard walls are used to transform the system from the periodic 
boundary conditions into the hard wall boundary conditions, as Davidchack and Laird did 
1^ . The detail is illustrated in Figure [TJ The two walls can be displaced synchronously 
along z direction from the boundary of the system (shown by the dot line in Figure [1]) to 
the position ai and a2 (|ai| = \a2 — Zq\ > cr/2, Zq is the dimension of the system in the 
2;-direction) , respectively. The interaction between the walls and the particles is 



where zl and are the coordinates of the particles and the walls in z direction, respectively. 
Therefore, when the separation between the walls and the boundaries is larger than the hard 
sphere radius, the walls do not interact with the particles and the system is in the bulk phase; 
however, when the walls are placed at the boundaries, the system are completely confined 
by the two hard walls. 




oo, when \zi — zZ,\ < <j/2, 
0, otherwise, 



(2) 
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The total Hamiltonian of the system if is a sum of the kinetic energy and the potential 
energy, given by 

^ = y — (3) 

i 

and 

(4) 

i<j 

where m is the mass of particles and Pi is the linear momenta. The Helmholtz free energy 
of the system is written as 

F{T, V, N, z^i, z^2) = -ksT ln(Q(T, V, N, z^i, 2^2)), (5) 

where Q(T, V, N, 

Zwi, Zu,2) is the canonical partition function of the system 

))• (6) 

In the leading order approximation, the total IFE is just the difference of the total free 



energy between the confined system and the bulk system 
is expressed by 

Fconfined{T, 7]) — Fbulk(T, T]) 
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23(1 . The IFE per unit area 7 



A 

F{T, 1], Zyji = 0, z^2 = Zq) - F{T, rj, z^i = al, z^2 = a2) 

A 



(7) 



here A is the total contact area between the system and two hard walls. 



Using WL algorithm [2l|, initially developed to compute the density of state(DOS) g{E), 
we can easily obtain the excess Helmholtz free energy in a Metropolis-type MC simula- 
tion. The original algorithm does a random walk in energy space with an acceptance ratio 
min{l, g{Ei)/ g{Ef)}, where Ei and Ej are respectively the energies of the current and the 
proposed configuration. The DOS of the accepted energy then adjusted by multiply it with 
a factor / > 1. The histogram is accumulated in the simulation and its flatness is tested. 
A new run with reduced / starts when the histogram reaches a specified flatness condition. 
The DOS obtained when the / reduced close to unity. 

The same procedure can be used to calculated the excess Helmholtz free energy. In 
order to use the method, we employ an expanded canonical ensemble, where the separation 
D^h from the hard walls to the system boundary is an additional ensemble variable. The 
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expanded ensemble method has been used to calculate the Helmholtz free energy at a series 



of temperatures [24| and the chemical potential of polymers [25| and colloid 26|]. Obviously, 
< -D^fe < Dq, where Dq = \ai\ = \a2 — Zq\. We use D^i, to rewrite ([7]) as 

1{V,T) = J . (8) 

Every D^b corresponds to a macroscopic state. The DOS g{Du,b) of the system can be 
computed via WL algorithm. For the hard interaction system, the logarithm of DOS is 
proportional to the partition function of the system. According to (j5]) and ([H]) we get 

^ kBT\n{g{D^b = Dp)) - kBTln{g{D^b = 0)) 

The expanded canonical ensemble simulations are performed in a similar way as in the 
canonical ensemble. In addition to the particles translation operation which is accepted or 
rejected in the usual way, the trial moves also include the walls displacement along z axis, 
which is accepted with the probability 

Pacc{Dt^ - DZr} = min{l, giD^^t) /giKin}, (10) 

Where and D^^l^ are the distances before and after the hard walls are displaced, respec- 
tively. In other words a random walk is implemented in D^b space instead of energy space. 
The calculation is as follows, the interval [0, -Do] is divided into many sub intervals, a random 
walk in this space is performed in the same way as that of the Wang-Landau method in the 
energy space. A initial guess of the DOS is assigned, the walk goes on with the acceptance 
probability given by (Q and a initial factor / > 1, the DOS of the accepted state is modified 
by multiply the corresponding DOS with /. A histogram of accepted states is accumulated 
and the flatness of it is tested. The run is finished when the histogram is fiat enough, here 
we use the criteria that the deviation of the histogram away from its average is less than 
20%. A new run starts with the DOS obtained in the previous run as initial DOS and the 
factor / is reduced to y/f. The full calculation is ended when the reduce factor / is very 
close to 1, here we use / — 1 < 10^^ — 10^^. At the end of the simulation the 7(77, T) can be 
determined from equation ([9]). 

III. RESULTS AND DISCUSSION 

The simulations were performed in a cuboid box with conventional periodic boundary 
conditions. The D^b is a continuous variable, however the WL algorithm demands a dis- 



Crete Duib- In the calculation the continuous is replaced with a set of discrete points 
and the walls may displace by jumping between nearest-neighbor points. Increasing D^b 
and decreasing Z)^;, are attempted with equal frequency, which are accepted or rejected by 
transition rule flTU]) . At first, we check the validity of our method by calculating and 
'jsy^ at two ordinary densities. 

A. Fluid-wall 7 

In this subsection we report the results of the 7 for the fluid-wall interface. The transverse 
dimensions of the simulation box are respectively = 8.1a and Ly = lo. In order to 
study the longitudinal flnite size effect we consider several different longitudinal dimensions 
Lz = 6.6(7, Lz = 13.2cr, Lz = 19. 80", = 26Aa and Lz = 33. la, and the number of particles 

are respectively 216, 432, 648, 864 and 1080, which correspond to a packing fraction of 
T] = 0.3. The dependence of 7^^ = 'jfwO''^ on Lz is shown in Figure El The 7^^ decreases 
with increasing the longitudinal dimension of the system at flxed packing fraction and the 
curve flattens gradually. In other words, there exists an effective interaction potential $eff 
between the two walls, which is a function of Lz and can be written as 

^os{Lz) = A{-ff^{Lz) - 7/««(oo)), (11) 

here A is the total contact area and 7^^(00) is the exact IFE. 

Our result of the 7^^ at the largest separation, Lz = 33.1, is 0.652 ± 0.005)kT/a^, which 
can be regarded as a good approximation of the exact IFE and this value is in very good 
agreement with that of Heni's [11], 7^^ = (0.656 ± 0.03)kT. We also performed simulations 
for a system with x ^ 2(Lx x Ly) and L'^ = 33. la, the results are basically the same 
as that of the results obtained from the smaller system which indicates our system size is 
already big enough. 

B. Solid- wall 7 

In the solid-wall case the initial conflguration is an ideal fee lattice with periodic boundary 
conditions. The system dimensions L^, Ly, and Lz are given in Tableland the number of 
particles A^ are determined in such a way that the volume fraction remains constant, rj = 0.6. 
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FIG. 2: Data points (with error bars) show the reduced fluid- wall interfacial free energy at 
r] = 0.3. The solid line is merely a guide to the eye. 
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TABLE I: The system dimensions and interfacial free energys with their statistical error for three 
crystal face orientations. 7*: our simulation result; 7J: MC simulation form Ref. [11]. 

Here, we calculate the crystal- wall 7s«, for the three crystalline orientations (1,1,1), (1,1,0) 
and (1, 0, 0). The results are compared with those of Heni's in Table [H The resulting 7*^ as 
a function of Lz are plotted in Figure [31 The profiles of the IFE and the effective potential 
are similar to those of the fluid phase. However, the IFEs are largely anisotropic and their 
ratio is 7110 : 7100 '■ 7iii ~ 5 : 3 : 2 at 77 = 0.6, which implies a hard sphere solid prefers to 
pack in the (1, 1, 1) orientation near a smooth hard wall. We also checked the lateral finite 
size effects of the crystal-wall '-jg^ for three orientations, all results are consistent with those 
obtained from the smaller systems within statistical error. 

C. Wetting 

The IFEs of the coexisting phases 7^^ and 7^^ are calculated in this section. Based on the 
results of the last subsections, we expect that the simulation systems with lateral dimensions 
La; X Ly ^ 60cr^ and longitudinal dimension Lz — 38o- is large enough to effectively reduce 
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FIG. 3: Data points show the crystal-wall reduced interfacial free energy 7*^ for three different 
orientations at = 0.6. From bottom to top: (1,1,1), (1,0,0) and (1,1,0) orientation. The solid 
line is merely a guide to the eye. Statistical errors do not exceed the symbol size. 

the finite-size effect. It is worth emphasizing that the coexisting volume fractions of the hard 
sphere fluid and solid phases are rjf = 0.4917 and rjs = 0.543 Q, [itI .^ 



rr 



27j, respectively, 



which are slightly different from those obtained by Hoover and Ree [2|. Because the hard 
sphere particles usually crystalize on a smooth hard wall with the (1, 1, 1) orientation, we 
need only consider this packing structure. The 7^^ can be easily determined, however, the 
prefreezing phenomenon prevents us from directly achieving the 7^^ from the simulation. 
We have to extrapolate the 7/„, from lower packing fraction to the bulk coexisting packing 
fraction rjf In order to get a more accurate estimation of the '-ff^ at the bulk coexisting 
packing fraction, we should use the data at packing fractions as close to the bulk coexisting 
as possible, on the other hand, the data used should correspond to fluid phases, i.e., before 
the prefreezing occurs. Therefore, it's crucial to accurately locate the volume fraction just 
before the prefreezing. Our simulation results are the relative Helmholtz free energy of 
the system with various Dy^b, therefore, if the prefreezing occurs so that crystalline layers 
appear at some D^b value, a singularity may be observed on the free energy versus D^h 
curve. Figure H] depicts the free energy-D^^ curve for different packing fractions rj close to 
the bulk coexisting packing. The volume fractions of the system from top to bottom are 
respectively T]f, 0.9975?7/, 0.9957]/, 0.9925?]/ and 0.99rif. From Figure H] we see that the 
packing fraction of forming the crystalline layers is lower than 0.995///, however it is difficult 
to locate the point exactly. With these observations it is appropriate to extrapolate from 
the volume fraction rj = 0.99?]/. On the other hand, by observing the density profile near the 




00 ' 02 ' 04 ' 06 ' 08 ' 10 



FIG. 4: (color online)The relative Helmholtz free energy versus the parameter D^f, for different 
packing fractions. From top to bottom the packing fraction rj = rjf, 0.9975r/j, 0.9957] f, 0.9925r]f 
and 0.99r]f. The shape of the upper curves obviously differs from the lower curves, which indicates 
the formation of the crystalline layers. The curves are shifted vertically for clarity. 

walls Dijkstra p] suggested that at the packing fraction r] ~ 0.9877]/ the crystalline layers 
begin to form. To check the consistence of the extrapolation, we also obtained the IFE by 
extrapolations from 'jjy^ at the volume fraction t] = O.dSdrjf = 0.4843. The two extrapolation 
results are the same within statistical error of the simulation. The results of the IFE are 
shown in Table HTl In Figure [5] we compare our simulation results with those from the scaled 



particle theory combined with the Carnahan-Starling equation of state(SPT) pjj. Here we 
don't present a comparison with Heni's simulation results, as the coexistence densities taken 
in their simulation are different from ours. 

Using the available latest crystal/melt IFE = 0.546(16) from simulation |l^ combined 
with our results, the contact angle 9 can be calculated in terms of the Young equation, cos 9 = 
il'fw~l'sw) ll% = 1.010(31), which nearly indicates that a complete wetting phase transition 
will occur. If the above is replaced by the value 0.55(2) obtained from nucleation 



experiment [28|, which is the average of the crystal/melt IFEs over all three orientations and 
larger than the 7^^ for (1,1,1) interface, we will have the same conclusion. The main source 
of error is the crystal/melt IFE 7^^,, a more accurate 7^^ is necessary to draw a definite 
conclusion. Of course, one can also improve the precision of the 7^^ and 7^^ by employing 
other sampling method, e.g. multicanonical algrithm 20 1. 
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0.0011 
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1.9105 


0.0016 


0.4867 


1.9323 


0.0018 


0.4917/ 


1.9872 


0.0012 


0.4917/1 


1.9896 


0.0016 


0.5430, 


1.4380 


0.0050 



TABLE II: The reduced surface free energys 7* with its error 6 for different volume fractions. The 
subscript / and fi refer to the surface tension of coexisting hard sphere fluid phase jj*^ obtained 
by extrapolating from 0.9977/ a-^icl 0.985///, respectively. And the subscript s refers to the surface 
free energy of coexisting hard sphere solid phase with (1, 1, 1) face. 

IV. CONCLUSION 

In conclusion, wc proposed a simple Monte Carlo method and calculated the IFEs of 
the hard sphere system at the smooth hard wall with high accuracy. The contact angle 
calculated with the accurate IFEs 7/^ and 7^^ make us confident that a smooth hard wall 
can be completely wetted by the hard sphere crystal at the wall-hard sphere fluid interface. 
Our conclusion is consistent with Dijkstra's work. However, when using the Young equation 
to study the wetting behavior, the finite-size effect can not be avoided completely as Dijkstra 
did in Ref. [6]. Therefore larger system and higher precision are necessary to confirm 
unambiguously the complete wetting phenomenon. Furthermore, the IFEs obtained also 
provide a benchmark to test other theoretical approaches. Our method can be extended 
to other systems (e.g. soft sphere system, polydisperse hard sphere system and hard/soft 
wall system) in a straightforward manner. We have already obtained some results of the 
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FIG. 5: The fluid-wall reduced surface tension versus the volume fraction. The solid line denotes 
SPT; the dashed line denotes a numerical fit from SPT; the filled square refers to our simulation 
results and * to the obtained from extrapolation. 



polydisperse hard sphere system which will be reported elsewhere. 

Note added in proof. — After the completion of this study, we are aware of Laird et al 
studied the same problem with molecular dynamics simulation and have the same conclusions 



as ours 
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